## This code creates Table 1, Figure 7 (a) and Figure 7 (c)

neighbor_dist_cat <- seq(10,100,10)

# Table 1, 1912 results

data <- colonial_towns_shp[colonial_towns_shp$neighbor_dist <= 50 & colonial_towns_shp$year==1912,]

coercive_specs <- c("coercive_OLS_BAGO_NoGEO","coercive_OLS_BAGO_NoFE","coercive_OLS_BAGO_FULL","coercive_OLS_BAGO_ih","coercive_OLS_BAGO_sd")

coercive_OLS_BAGO_NoGEO <- conleyreg(coercive~replaced+log(houses),data=data@data,dist_cutoff=100,lat="latitude",lon="longitude",gof=TRUE)
coercive_OLS_BAGO_NoFE <- conleyreg(coercive~replaced+log(houses)+altitude+slope_idx+as.factor(soil_quality)+log(yangon_dist)+log(river_dist)+longitude_2+latitude_2,data=data@data,dist_cutoff=100,lat="latitude",lon="longitude",gof=TRUE)
coercive_OLS_BAGO_FULL <- conleyreg(coercive~replaced+log(houses)+altitude+slope_idx+as.factor(soil_quality)+log(yangon_dist)+log(river_dist)+as.factor(district)+longitude_2+latitude_2,data=data@data,dist_cutoff=100,lat="latitude",lon="longitude",gof=TRUE)
coercive_OLS_BAGO_ih <- conleyreg(coercive_ih~replaced+log(houses)+altitude+slope_idx+as.factor(soil_quality)+log(yangon_dist)+log(river_dist)+as.factor(district)+longitude_2+latitude_2,data=data@data,dist_cutoff=100,lat="latitude",lon="longitude",gof=TRUE)
coercive_OLS_BAGO_sd <- conleyreg(coercive_sd~replaced+log(houses)+altitude+slope_idx+as.factor(soil_quality)+log(yangon_dist)+log(river_dist)+as.factor(district)+longitude_2+latitude_2,data=data@data,dist_cutoff=100,lat="latitude",lon="longitude",gof=TRUE)

coercive_table <- matrix(nrow = 7,ncol=6)

for(i in 1:length(coercive_specs)){
  coercive_table[1,i+1] <- coercive_specs[i]
  model <- get(coercive_specs[i])
  coef <- round(model$coefficients[2],digits=4)
  se <- round(model$coefficients[2,2],digits=4)
  coercive_table[2,i+1] <- coef
  coercive_table[3,i+1] <- paste0("(",se,")")
  coercive_table[4,i+1] <- paste0("[",round(coef-1.96*se,digits=4),",")
  coercive_table[5,i+1] <- paste0(round(coef+1.96*se,digits=4),"]")
  coercive_table[6,i+1] <- nobs(model)
  coercive_table[7,i+1] <- round(model$adj.r.squared,digits=4)
}

xtable(coercive_table)

# Figure 7 (a)

coercive_BAGO <- data.frame(coefficient=double(),se=double(),upper=double(),lower=double())

for(j in 1:length(neighbor_dist_cat)){
  model <- conleyreg(coercive_sd~replaced+log(houses)+altitude+slope_idx+as.factor(soil_quality)+log(yangon_dist)+log(river_dist)+as.factor(district)+longitude_2+latitude_2,dist_cutoff=100,lat="latitude",lon="longitude",
                     data=colonial_towns_shp@data[colonial_towns_shp@data$year==1912 & colonial_towns_shp@data$neighbor_dist<=neighbor_dist_cat[j],])
  coercive_BAGO[j,"coefficient"] <- model[2,1]
  coercive_BAGO[j,"se"] <- model[2,2]
  coercive_BAGO[j,"upper"] <- coercive_BAGO[j,"coefficient"] + (1.96*coercive_BAGO[j,"se"])
  coercive_BAGO[j,"lower"] <- coercive_BAGO[j,"coefficient"] - (1.96*coercive_BAGO[j,"se"])
  coercive_BAGO[j,"upper_10pc"] <- coercive_BAGO[j,"coefficient"] + (1.645*coercive_BAGO[j,"se"])
  coercive_BAGO[j,"lower_10pc"] <- coercive_BAGO[j,"coefficient"] - (1.645*coercive_BAGO[j,"se"])
  coercive_BAGO[j,"pch"] <- ifelse(sign(coercive_BAGO[j,"upper"])==sign(coercive_BAGO[j,"lower"]),"black",ifelse(sign(coercive_BAGO[j,"upper_10pc"])==sign(coercive_BAGO[j,"lower_10pc"]),"gray","white"))
}

x.axis <- rep(1:10)
labels <- seq(0,100,10)
labels.rev <- rev(labels)
png(file=paste0("Plots/coercive_BAGO_1912.png",sep=""), width=6.5, height=6.5, pointsize=9, units="in", res=600)
par(mar=c(6, 6, 3, 3) + 0.1)
plot(c(0.5,10.5),c(-0.5,0.5),type="n",axes=F,
     xlim=range(0.5:10.5),
     ylim=range(-0.5:0.5),
     xlab="Distance from a location in sittan (km)",
     ylab="Coefficient of headman replacement",cex.lab=1.5,cex.axis=1.5)
box()
axis(2, at=seq(-0.5,0.5,by=0.1), labels=c(-1,-0.8,-0.6,-0.4,-0.2,0,0.2,0.4,0.6,0.8,1))
axis(1, at=seq(0,10), las=2, labels=labels)
segments(x.axis,coercive_BAGO$lower/2,x.axis,coercive_BAGO$upper/2)
points(x.axis,coercive_BAGO$coefficient/2, pch=21, bg=coercive_BAGO$pch)
abline(h=0, lty=3)
dev.off()

# Table 1, 1924 results

data <- colonial_towns_shp[colonial_towns_shp$neighbor_dist <= 50 & colonial_towns_shp$year==1924,]

data$longitude_2 <- data$longitude

data$latitude_2 <- data$latitude

coercive_specs <- c("coercive_OLS_BAGO_NoGEO","coercive_OLS_BAGO_NoFE","coercive_OLS_BAGO_FULL","coercive_OLS_BAGO_ih","coercive_OLS_BAGO_sd")

coercive_OLS_BAGO_NoGEO <- conleyreg(coercive~replaced+log(houses),data=data@data,dist_cutoff=100,lat="latitude",lon="longitude",gof=TRUE)
coercive_OLS_BAGO_NoFE <- conleyreg(coercive~replaced+log(houses)+altitude+slope_idx+as.factor(soil_quality)+log(yangon_dist)+log(river_dist)+longitude_2+latitude_2,data=data@data,dist_cutoff=100,lat="latitude",lon="longitude",gof=TRUE)
coercive_OLS_BAGO_FULL <- conleyreg(coercive~replaced+log(houses)+altitude+slope_idx+as.factor(soil_quality)+log(yangon_dist)+log(river_dist)+as.factor(district)+longitude_2+latitude_2,data=data@data,dist_cutoff=100,lat="latitude",lon="longitude",gof=TRUE)
coercive_OLS_BAGO_ih <- conleyreg(coercive_ih~replaced+log(houses)+altitude+slope_idx+as.factor(soil_quality)+log(yangon_dist)+log(river_dist)+as.factor(district)+longitude_2+latitude_2,data=data@data,dist_cutoff=100,lat="latitude",lon="longitude",gof=TRUE)
coercive_OLS_BAGO_sd <- conleyreg(coercive_sd~replaced+log(houses)+altitude+slope_idx+as.factor(soil_quality)+log(yangon_dist)+log(river_dist)+as.factor(district)+longitude_2+latitude_2,data=data@data,dist_cutoff=100,lat="latitude",lon="longitude",gof=TRUE)

coercive_table <- matrix(nrow = 9,ncol=6)

for(i in 1:length(coercive_specs)){
  coercive_table[1,i+1] <- coercive_specs[i]
  model <- get(coercive_specs[i])
  coef <- round(model$coefficients[2],digits=4)
  se <- round(model$coefficients[2,2],digits=4)
  coercive_table[2,i+1] <- coef
  coercive_table[3,i+1] <- paste0("(",se,")")
  coercive_table[4,i+1] <- paste0("[",round(coef-1.96*se,digits=4),",")
  coercive_table[5,i+1] <- paste0(round(coef+1.96*se,digits=4),"]")
  coercive_table[6,i+1] <- nobs(model)
  coercive_table[7,i+1] <- round(model$adj.r.squared,digits=4)
}

xtable(coercive_table)

# Figure 7 (c)

coercive_BAGO <- data.frame(coefficient=double(),se=double(),upper=double(),lower=double())

for(j in 1:length(neighbor_dist_cat)){
  model <- conleyreg(coercive_sd~replaced+log(houses)+altitude+slope_idx+as.factor(soil_quality)+log(yangon_dist)+log(river_dist)+as.factor(district)+longitude_2+latitude_2,dist_cutoff=100,lat="latitude",lon="longitude",
                     data=colonial_towns_shp@data[colonial_towns_shp@data$year==1924 & colonial_towns_shp@data$neighbor_dist<=neighbor_dist_cat[j],])
  coercive_BAGO[j,"coefficient"] <- model[2,1]
  coercive_BAGO[j,"se"] <- model[2,2]
  coercive_BAGO[j,"upper"] <- coercive_BAGO[j,"coefficient"] + (1.96*coercive_BAGO[j,"se"])
  coercive_BAGO[j,"lower"] <- coercive_BAGO[j,"coefficient"] - (1.96*coercive_BAGO[j,"se"])
  coercive_BAGO[j,"upper_10pc"] <- coercive_BAGO[j,"coefficient"] + (1.645*coercive_BAGO[j,"se"])
  coercive_BAGO[j,"lower_10pc"] <- coercive_BAGO[j,"coefficient"] - (1.645*coercive_BAGO[j,"se"])
  coercive_BAGO[j,"pch"] <- ifelse(sign(coercive_BAGO[j,"upper"])==sign(coercive_BAGO[j,"lower"]),"black",ifelse(sign(coercive_BAGO[j,"upper_10pc"])==sign(coercive_BAGO[j,"lower_10pc"]),"gray","white"))
}

x.axis <- rep(1:10)
labels <- seq(0,100,10)
labels.rev <- rev(labels)
png(file=paste0("Plots/coercive_BAGO_1924.png",sep=""), width=6.5, height=6.5, pointsize=9, units="in", res=600)
par(mar=c(6, 6, 3, 3) + 0.1)
plot(c(0.5,10.5),c(-0.5,0.5),type="n",axes=F,
     xlim=range(0.5:10.5),
     ylim=range(-0.5:0.5),
     xlab="Distance from a location in sittan (km)",
     ylab="Coefficient of headman replacement",cex.lab=1.5,cex.axis=1.5)
box()
axis(2, at=seq(-0.5,0.5,by=0.1), labels=c(-1,-0.8,-0.6,-0.4,-0.2,0,0.2,0.4,0.6,0.8,1))
axis(1, at=seq(0,10), las=2, labels=labels)
segments(x.axis,coercive_BAGO$lower/2,x.axis,coercive_BAGO$upper/2)
points(x.axis,coercive_BAGO$coefficient/2, pch=21, bg=coercive_BAGO$pch)
abline(h=0, lty=3)
dev.off()

